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MANPAKs  A  SET  OF  ALGORITHMS  FOR 
COMPUTATIONS  ON  IMPLICITLY  DEFINED  MANIFOLDS1 


BY 

Werner  C.  Rheinboldt* 

Abstract.  Mathematical  models  often  involve  differentiable  manifolds  that  are  implicitly 
defined  as  the  solution  sets  of  systems  of  nonlinear  equations.  The  resulting  computational 
tasks  differ  considerably  from  those  arising  for  manifolds  defined  in  parametric  form.  Here 
a  collection  of  algorithms  is  presented  for  performing  a  range  of  essential  tasks  on  general, 
implicitly  specified  submanifolds  of  a  finite  dimensional  space.  This  includes  algorithms  for 
determining  local  parametrizations  and  their  derivatives,  and  for  evaluating  quantities  related 
to  the  curvature  and  with  sensitivity  measures.  The  methods  have  been  implemented  as  a 
FORTRAN  77  package,  called  MANPAK. 


1.  Introduction. 

Mathematical  models  of  many,  practically  important  scientific  and  technical  problems 
involve  differentiable  manifolds  that  are  implicitly  defined  as  the  solution  sets  of  systems 
of  nonlinear  equations.  For  example,  the  computational  study  of  equilibria  typically  leads 
to  nonlinear  systems  of  the  form 

(U)  F(z.  A)  =  0.  F:SmxR^r, 

where  F  is  a  sufficiently  smooth  mapping,  r  <=  Rm  a  state  variable,  and  A  €  R*  a  parameter 
vector.  Here,  interest  often  centers  on  d.  termining  the  behavior  of  the  solutions  under 
variation  of  A.  Under  simple  conditions  (see  Section  2  below),  the  zero  set  M  =  {(*,A)  6 
Rm  x  R'*  :  F(z.  A)  =  0}  has  the  structure  of  a  submanifold  of  dimension  d  of  the  product 
Rm  x  R*  of  the  state  and  parameter  space.  Then  we  are  faced  with  a  computational  analysis 
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of  particular  features  of  this  manifold  XI,  such  as.  for  instance,  of  certain  types  of  singular 
points  on  M  or  of  the  curvature  behavior  of  XI.  Another  example  arises  in  connection 
with  equality  constrained  dynamical  systems  that  are  modelled  by  differential-  algebraic 
equations  (DAEs).  Such  DAEs  are  known  to  be  closely  related  to  ordinary  differential 
equations  (ODEs)  on  implicitly  defined  differentiable  manifolds. 

For  these,  and  similar  problems,  efficient  numerical  methods  are  required  for  computa¬ 
tions  on  implicitly  defined  submanifolds  of  Rn.  These  tasks  differ  considerably  from  those 
encountered,  for  example,  in  computer  graphics,  where  curves  or  surfaces,  i.e.  submani¬ 
folds  of  R3,  are  considered  that  are  explicitly  specified  in  a  parametric  form  M  =  {x  € 
R3  :  x  =  v>(y),  y  €  R**}  with  d  =  1  or  d  =  2,  respectively.  In  fact,  one  of  the  basic  compu¬ 
tational  problems  arising  in  connection  with  any  implicitly  defined  manifold  is  exactly  the 
construction  of  such  parametrizations  and  their  derivatives  which  requires  the  solution  of 
certain  systems  of  nonlinear  equations.  Other  related  problems  concern,  for  instance,  the 
computation  of  the  curvature  of  the  manifolds  or  of  the  sensitivity  of  the  solutions  under 
specific  changes. 

This  paper  presents  numerical  methods  for  performing  these  and  related  tasks  on  gen¬ 
eral,  implicitly  specified  submanifolds  of  a  finite  dimensional  space.  The  methods  have 
been  collected  in  a  FORTRAN  77  package,  called  MANPAK.  All  routines  use  reverse  com¬ 
munication  to  avoid  calls  to  subroutines  for  the  evaluation  of  user-defined  functions.  The 
package  is  intended  for  applications  to  small  or  medium-sized  problems,  mainly  because 
they  involve  many  dense  matrix  computations.  Some  examples  for  use  of  the  algorithms 
are  noted  here,  and,  in  addition,  we  refer  to  the  companion  paper  [Rh96cj  for  applications 
to  differential-algebraic  equations. 

2.  Background. 

For  ease  of  reference,  this  section  collects  some  basic  definitions  and  theorems  about 
submanifolds  of  Rn.  For  details  and  proofs  see,  e.g.,  the  texts  [S79],  (AMRS8). 

Let  F  :  Rn  *-»  Rm  be  of  class  Ck ,  k  >  0,  on  an  open  set  E  C  R";  that  is,  assume  that  F 
is  continuous  and  that  ail  its  partial  derivatives  of  order  at  most  k  exist  and  are  continuous 
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on  £.  For  k  >  1.  £  is  an  immersion  or  submersion  at  a  point  x  €  E  if  its  first  derivathe 
DF(x)  6  £(Rn,Rm)  is  a  one-to-one  (linear)  mapping  or  a  mapping  onto  Rm.  respectively. 
We  call  F  a  submersion,  or  immersion  on  a  subset  5  of  £  if  it  has  that  property  at  each 
point  of  5.  These  definitions  obviously  require  that  n  <  m  for  F  to  be  an  immersion  and 
n  >  m  for  it  to  be  a  submersion.  Clearly,  if  n  >  m  and  DF(x)  has  maximal  rank  m.  then 

F  is  a  submersion  at  x. 

A  nonempty  subset  M  C  R"  is  a  submanifold  of  Rn  of  dimension  d  and  class  Ck  if 
for  every  xc  €  M  there  exists  an  open  neighborhood  Vn  of  xc  in  R"  and  a  submersion 
F  .  V"  ~  of  class  Ck  such  that  M  D  V"  =  {x  €  V"  :  F(x)  =  0}.  In  particular,  if 

F  ;  £  Rm,  n  -  m  =  d  >  0,  is  of  class  Ck  on  an  open  set  E  C  R"  and  a  submersion 
on  :=  {x  €  £  :  F{x)  =  0},  then  M  is  a  d-dimensional  Ck  submanifold  of  Rn.  Any 
non-empty,  (relatively)  open  subset  of  a  d-dimensional  C*- submanifold  of  Rn  is  itself  a  Ck 
submanifold  of  Rn  of  the  same  dimension. 

For  the  analysis  of  submanifolds  of  R"  we  need  local  parametrizations.  Let  M  be  any 
nonempty  subset  of  Rn.  A  local  d-dimensional  Ck  parametrization  of  M  is  a  pair 
consisting  of  a  nonempty,  open  subset  Vd  of  Rd  and  a  Cfc  mapping  '•  V*  ~  R"  such  that 
^(V*)  is  (relatively)  open  in  M,  *  is  a  homeomorphism  of  Vd  onto  it  image  v’lV'1),  and 
^  is  an  immersion  on  Vd.  For  any  point  ze  of  M  such  that  ic  €  <s(Vd)  we  call  ( Vd ,s?)  a 
local  d- dimensional  Cfc  parametrization  of  M  near  xc.  A  nonempty  subset  M  C  Rn  is  a 
d-dimensional  Ck  submanifold  of  Rn  if  and  only  if  for  every  xc  €  M  there  exists  a  local 
d-dimensional  Ck  parametrization  of  M  near  xc.  If  M  is  a  d-dimensional  Ck  submanifold 
of  Rn,  and  ( Vr,  \p)  a  local  r-dimensional  Ck  parametrization  of  A/,  then,  necessarily,  r  -  d. 

Instead  of  defining  tangent  spaces  in  general,  we  use  here  the  following  characterization: 
Let  M  be  a  d-dimensional  C*  submanifold  of  R».  For  any  point  xc  6  Af  we  can  choose,  by 
definition,  an  open  neighborhood  V"  C  Rn  of  xc  and  a  submersion  £  :  V"  -  R"-'  at  xc 
such  that  M  H  V"  =  {x  €  V"  :  F(x)  =  0}.  Then,  it  can  be  shown  that  the  d- dimensional 
linear  subspace  5  =  ker££(xc)  of  R"  is  independent  of  the  particular  choice  of  the  local 
submersion  £;  that  is,  5  depends  only  on  M  and  the  particular  point.  This  space  5  is  the 
tangent  space  of  M  at  xc  and  is  denoted  by  TItM.  The  subset  TM  = 
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of  Rn  x  Rn  is  the  tangent  bundle  of  M. 

Since  TtRn  =  Rn  for  every  x  €  R",  we  have  TRn  =  R"  x  R".  Thus,  the  tangent  bundle 
of  a  submanifold  M  of  Rn  appears  as  a  subset  of  TRn.  In  general,  TM  is  a  submanifold 
of  TRn  =  Rn  x  Rn.  More  specifically,  if  M  be  a  d-dimensional  Ck  submanifold  of  Rn  with 
k  >  2,  then,  TM  is  a  2d-dimensional  Ck~l  submanifold  of  TRn  =  R"  x  R"  2  RJn.  For 
k  >  2  the  local  parametrizations, of  the  Ck~l  submanifold  TM  of  R2n  can  be  constructed 
easily  from  local  Ck  parametrizations  of  M.  Let  M  be  a  d-dimensional  Ck  submanifold  of 
Rn,  k  >  2,  and  (xc,ve)  €  TM.  Then  for  any  local  Ck  parametrization  (V*,p)  of  M  near 
xcy  the  pair  (V*  x  R*,  (p,  D<p))  is  a  Cfc_1  local  parametrization  of  TM  near  (xc,ve). 

3.  Computation  of  Local  Parametrizations. 

This  section  presents  algorithms  for  computing  local  parametrizations  on  submanifolds 
of  Rn  which  are  implicitly  defined;  by  local  submersions.  Some  of  the  material  was  given 
earlier,  in  part,  in  (RhS8],  [Rh80],  [Rh96b].  In  view  of  the  local  nature  of  the  methods, 
there  is  no  loss  of  generality  to  restrict  attention  to  the  case  of  a  single  submersion. 

Assumption  A:  With  positive  integers  d,  k,  m,  n  such  that  m  =  n  —  d,  n  >  d,  let 
F  :  E  h*  Rm  be  a  Ck  mapping  on  an  open  subset  E  of  Rn  and  a  submersion  on  M  = 
jr-i(o)  =  [x  $  E  :  F(x)  =  0}  whence  M  is  a  d-dimensional  Ck  submanifold  of  Rn. 

~h  V-.  ;  i  . 

The  following  result  exhibits  a  method  for  the  computation  of  a  local  parametrization 

on  M : 

Theorem  3.1:  Under  Assumption  .4,  let  U  :€  £(R^,R")  be  a  linear  isomorphism  from 
R<*  onto  a  d-dimensional  linear  subspace  T  C  Rn.  Denote  by  U*  6  £(R",  R*)  the  adjoint 
of  U  and  by  J  :  R*  *-»  Rm  x  R*  the  canonical  injection  that  maps  R*  isomorphically  onto 
{0}  x  R*.  Then  the  Ck  mapping  H  :,E~  Rm  x  Rd,  defined  by  H(x)  =  (. F{x)yU*x )  for 
x  €  E.  is  a  local  diffeomorphism  near  a  point  xe  €  M  if  and  only  if 

(3.1)  T,tMnTL  =  {  0}. 
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If  (3.1)  holds  al  x.  then  there  exists  an  open  set  Vd  of  Rd  such  that  the  pair  (V  ■  r- ). 
defined  with  the  mapping  9  =  ff-  o  J  :  V“  ~  R\  is  a  iocai  C‘  parara.tmation  of  M 

near  xc. 

Proo/:  By  Assumption  A,  ff  is  of  class  Ck  on  E.  Evidently,  DH(xe)h  =  0  for  any  h  €  R" 
requires  that  DF{xc)h  =  0  and  =  0  whence.  A  €  TrcA/  and,  because  of 

(3.2)  (h,Uy)  «  (Umh,y)  **  0,  y  €  R*4, 

that  h  €  Tx  which  by  (3.1)  implies  that  h  =  0.  Conversely,  if  there  exists  a  nonzero 
/i  €  TteM  n  T1  then  DE(xe)h  -  0  and  (3.2)  requires  that  E'h  =  0  which  together 
shows  that  DH(xe)h  =  0.  Hence  there  is  some  open  neighborhood  Vn  of  xe  in  Rn  such 
that  if  is  a  diffeomorphisra  from  Vn  onto  the  open  set  H(Vn)  in  Rn.  Evidently,  the  set 
H(M  (lVn)  =  if(Vn)n({0}  x  R*)  isopen  in  {0}  x  R*  and  J~XH(MC\  Vn)  =  Vd  C  Rn  is 
an  open  subset  of  R-.  This  shows  that  9  »  H~x  o  J  is  a  C*  mapping  from  V*  onto  the 
open  subset  M  n  Vn  of  M.  Both  p  and  its  inverse  J“l  o  ff\Mnv*  are  continuous  and  hence 
*  is  a  homeomorphism  of  V*  onto  M  fi  V".  Since  both  H~'  and  J  are  immersions,  the 
same  is  true  for  p.  Thus,  altogether,  p  is  a  local  d-dimensional  Ck  parametrization  of  M  ^ 

near  xe.  □  '|| 

Note  that  (3.1)  is  equivalent  with  ker  DF(xc)  n  TL  =  {0}  or  rge  DF{xe)  n  T  =  {0}.  J| 
At  any  point  xc  €  M  an  obvious  choice  for  a  subspace  T  satisfying  (3.1)  is,  of  course, 
j  -  TttM.  For  ease  of  reference  we  introduce  the  following  terminology: 

Definition  3.3:  Under  Assumption  A  a  d-dimensional  linear  subspace  T  C  Rn  is  a  co- 
ordinate  subspace  of  M  at  xc  €  M  if  (3.1)  holds,  otherwise  x,  is  a  fcldpomt  of  M  with 
respect  to  T.  In  the  case  T  =  T„M  we  speak  of  the  tangential  coordinate  space  of  M  at 

the  point  xc. 

Theorem  3.1  readily  becomes  a  computational  procedure  for  local  parametnzations  by 
the  introduction  of  bases.  On  Rn  and  Rd  the  canonical  bases  e?, ...  ,eS  and  ef,.. .  ,en, 
respectively,  will  be  used  and  we  assume  that  the  vectors  u, - €  Rn  form  an  orthonor¬ 

mal  basis  of  the  given  coordinate  subspace  T  of  M  at  xe.  Then  the  matrix  representation 


5 


of  the  mapping  V  is  the  n  x  d  matrix  with  the  vectors  Uj, . . .  .  u<i  as  columns.  We  denote 
this  matrix  again  by  U.  It  is  advantageous  to  shift  the  open  set  Vd  such  that  <,3(0)  =  xe. 
Then,  in  component  form,  the  nonlinear  mapping  H  assumes  the  form 

(3.2)  H  :Rn  >-*  Rn,H(x)  =  (  ).  Vi  €  E  C  Rn 

\UT{x-xc)J 

where  F(x)  is  the  column  vector  consisting  of  the  m  components  of  F  evaluated  at  x.  By 
definition  of  9  we  have 


(3.3) 


H(v(y))  =  Jy,  Vy€Vd. 


Thus,  the  evaluation  of  x  =  9(y)  for  given  y  €  Vd  requires  the  solution  of  the  nonlinear 
system  of  equations 


(3.4) 


H(x)  = 


(  F(x) 
\UT(x-xe) 


Since  (3.1)  is  assumed  to  hold  at  xc  6  M,  the  Jacobian 


(3.5) 


DH(x) 


(T) 


is  nonsingular  in  an  open  neighborhood  of  x  —  xc. 

Experience  has  shown  that,  for  the  solution  of  the  nonlinear  system  (3.4),  a  chord 
Newton  method  works  well  in  practice.  It  is  advantageous  to  start  with  x°  =  xc  +  U y 
which  allows  the  process  to  be  applied  in  the  y-independent  form 

(F(x)  \ 

Q  J  ’  A  =  DH{ze). 

We  sketch  briefly  its  convergence  properties.  For  given  e  >  0  such  that  ||>1  Ml*  — 
there  exists,  by  Assumption  A,  a  6  >  0  such  that  the  closed  ball  B  s  B(xc6)  is  contained 
in  E,  and  that  ||£F(x)  -  £>F(xc)||2  <  *  for  all  x  €  B.  Then 

||^(i)ll2<P",||2  \\DF(xc)-DF(x)\\2  < 


Vx  €  B 


'hows  that  G  is  contractive  on  D.  Moreover,  for  x  €  B  it  follows  from 

||G(x)  -  G(xr)j|,  <  |iA_l||.,  ||  /  DFitrTMx -xM)(x-xe)ds|| 

Jo 

<  2  «||.4"'!Ux  -  jv!|,  1  lix-  ^||2 

that  G  maps  D  into  itself.  Recall  from  the  proof  of  Theorem  3.1  that  //  is  a  diffcomorphism 
from  an  open  neighborhood  V"  of  xr  onto  its  image.  Let  h  be  sufficiently  small  that 
B  C  Vn.  Then  it  follows  from  the  contraction  theorem  that,  for  any  y  €  Vd.  the  process 
(3.6),  started  from  x°  =  xc  +  Uy,  converges  to  the  unique  fixed  point  x*  €  D  of  C.  Clearly, 

F(x*)  =  0  and.  from 

UT(X  -  Xc)  =  UT(x°  -  xe)  +  f;  UT{XJ  +  1  -  X1)  =  y 

;=0 

we  obtain  that  H(x')  =  0  and  therefore,  in  view  of  x*  €  V'n.  that  x*  =  y>(y). 

This  shows  that,  for  any  local  vector  y  near  the  origin  of  the  following  algorithm 
produces  the  point  x  =  <p{y)  in  the  local  parametrization  (Vd,s)  near  xc  defined  by 

Theorem  3.1: 

GPHI  Input:  Center  point  x,.  €  M  <>f  the  local  parametrization,  local  vector  y  €  R  , 
basis  matrix  L  .  Jacobian  DF{xc)i  tolerances, 
x  :=  xr  +  V  <J\ 

Compute  the  LU  factorization  of  A  - 
While:  iterates  do  not  meet  tolerances 
Return  for  the  evaluation  of  F(x)\ 


Solve  Aw  =  q  for  w  €  Rn: 
x  :=  x  —  xv. 

End  While 
Output:  y-(y)  :=  x- 

For  the  sake  of  clarity,  our  reverse  communication  paradigm  was  here  only  indicated 
bv  the  statement  ‘Return  for  the  evaluation  of  F(xf  The  MANPAK  implementation  of 


DF(xe)\ 

v  r 


GPHI  uses  a  different  step-order  for  handling  the  repeated  returns  to  the  calling  program. 
Of  course,  a  full  Newton  method  or  some  other  iterative  process  can  be  applied  as  well. 
However,  experience  has  shown  that  the  faster  convergence  of  Newton’s  method  at  the 
expense  of  several  evaluations  and  factorizations  of  DH  does  not  improve  the  overall 
performance. 

So  far,  it  was  assumed  that  T  is  a  given  coordinate  subspace  of  M  at  xc  and  that  an 
orthonormal  basis  of  T  is  available.  For  the  construction  of  T  we  use  a  simple  reformu¬ 
lation  of  (3.1):  Let  Z  =  TL  be  the  orthogonal  complement  of  T  (under  the  canonical 
inner  product  of  Rn),  then  (3.1)  implies  that  Z  ©  Tt<M  -  Rn;  that  is,  that  Z  is  also 
a  complementary  subspace  of  TteM.  Thus,  in  order  to  ensure  the  validity  of  (3.1),  it  is 
advantageous  to  construct  T  by  choosing  a  complementary  linear  subspace  Z  of  Tt<  M  and 
to  determine  T  as  the  orthogonal  complement  of  Z . 

As  before  we  use  the  canonical  bases  on  R^  and  Rn  and  assume,  for  the  moment,  that 
a  basis  zx , . . .  ,  zm  of  Z  is  available.  Then,  T  is  the  nullspace  of  the  n  x  m  matrix  formed 
with  these  vectors  as  its  columns.  We  denote  this  matrix  again  by  Z.  Our  task  is  now 
to  compute  an  orthogonal  basis  of  the  nullspace  of  the  transposed  matrix  ZT  for  which 
obviously  rank  ZT  =  m.  There  are  several  approaches  for  this;  probably  the  simplest  one 
is  based  on  the  LQ-factorization  (with  row  pivoting) 

(3.7)  Zt  =  Pt(L  0  )Qr,  0  =  (<?i,<?2). 

Here  P  is  an  m  x  m  permutation  matrix,  L  an  m  x  m,  nonsingular,  lower  triangular  matnx, 
and  Qannxn  orthogonal  matrix  partitioned  such  that  Q\  and  Q2  are  n  x  m  and  n  x  d 
matrices,  respectively.  Then,  clearly,  the  d  columns  of  Q2  form  the  desired  orthonormal 
basis  of  T.  This  justifies  the  following  MANPAK  algorithm: 

COBAS  Input:  m  x  n  matrix  ZT  of  rank  m; 

Compute  the  LQ-factorization  (3.7)  of  ZT  with  row-pivoting; 

For  j  =  1 . d  Do:  u;  :=  Qjef, 

Output:  U  ,u*). 
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Other  algorithms  for  the  computation  of  nuilspace-bases  of  m  x  n  matrices  are  given, 
for  example,  in  [BHKS5],  [CPS6]. 

Obviously,  when  the  tangential  coordinate  system  is  used  at  xe,  then  COBAS  can  be 
applied  with  the  Jacobian  matrix  DF(xc)  as  the  matrix  ZT .  In  that  case,  GPHI  simplifies 
considerably  if  the  LQ-factorization  (3.7)  of  DF(xr)  is  applied  for  the  solution  of  the 
corrector  equation  Aw  =  q  in  GPHI.  In  fact,  this  equation  has  the  block-components 
DF(x)w  =  F(x)  and  UTw  =  0,  which,  with  (3.7)  and  U  =  Q2,  can  be  rewritten  as 
LUtw  =  PF{x)  and  UTw  =  0.  Thus,  in  this  case,  the  algorithm  can  be  modified  as 

follows: 

TPHI  Input:  Center  point  re  €  M  <>f  the  local  parametrization.  local  vector  y  €  R*. 
the  LQ  factorization  (3.6)  of  DF(xc),  tolerances: 
x  :=  xc  +  Qj y; 

While:  iterates  do  not  meet  tolerances 
Return  for  the  evaluation  of  F(x)\ 

Solve  Lv  =  PF(x)  for  v  €  Rm; 
x  :=  x  -  Qjv, 

End  While 

Output:  v?(y)  :=  x. 

Thus,  for  each  iteration  step  we  need  to  solve  now  only  an  m  x  m,  rather  than  an  n  x  n, 
lower-triangular.  The  convergence  behavior,  of  course,  remains  the  same. 

In  COBAS  the  transposed  basis  matrix  Zr  of  some  complementary  space  of  TtcM  was 
assumed  to  be  available.  In  order  to  compute  such  a  matrix  it  is  natural  to  start  with  the 
matrix  respresentation  DF(xc)  of  the  Jacobian  of  F  at  xc  for  which  the  rows  form  a  basis 
of  {TZeM)L-  We  construct  bases  of  complementary  spaces  of  TteM  by  replacing  suitable 

columns  vectors  of  DF{xc). 

In  applications  it  is  frequently  important  to  work  with  coordinate  spaces  T  that  contain 
a  specific  canonical  basis  vector,  say.  e?,  of  Rn.  Often  the  reason  for  this  is  that  the 
independent  variable  rt  represented  by  this  vector  is  of  a  special  nature,  as,  for  instance, 
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when  X(  corresponds  to  the  time  variable  in  some  non-autonomous  DAE.  Evidently,  in 
order  to  ensure  that  t\ '  €  T  we  have  to  replace  the  £-th  column  of  DF(xe )  by  a  zero 
column.  This  leads  to  the  following  MANPAK  algorithm: 

GNBAS  Input:  The  Jacobian  matrix  DF(xc),  index  £  €  {1,...  ,n}; 

Form  Zr  by  zeroing  column  t  of  DF(xc) ; 

Use  COBAS  to  compute  the  orthonormal  basis  U  of  of  ker  W\ 

OUTPUT  The  basis  matrix  U. 

Of  course,  this  algorithm  requires  that  the  constructed  matrix  ZT  still  has  maximal 
rank  m.  In  order  to  verify  this,  while  computing  the  basis  of  the  nullspace,  we  may  replace 
the  LQ-factorization  of  ZT  in  COBAS  by  the  singular  value  decomposition  (SVD)  and 
then  apply  scaling  and  standard  rank-tests  (see  e.g.  [GoVL89]).  Such  a  version  of  the 
algorithm  was  not  included  in  MANPAK  since  it  appears  to  be  rarely  needed. 

In  certain  applications  it  is  desirable  to  work  with  coordinate  spaces  spanned  by  d 
suitably  chosen  canonical  basis  vectors  e* , . . .  , e",  of  Rn,  (see  e.g.,  [WH82],  [PoRh91]).  In 
this  case  the  space  ZT  should  be  spanned  by  n  -  d  canonical  basis  vectors  of  Rn.  These 
vectors  can  be  obtained  from  the  QR-factorization  of  the  Jacobian  DF(xc)  with  column 
pivoting.  The  resulting  permutation  selects  m  column- vectors  of  the  Jacobian  such  that 
the  m  x  m  matrix  formed  by  these  columns  is  nonsingular.  The  indices  of  these  selected 
columns  correspond  to  the  desired  m  canonical  basis  vectors  of  Rn  spanning  the  space  Z. 
This  leads  to  the  following  algorithm: 

GCBAS  Input:  The  Jacobian  matrix  DF(x c) 

Compute  the  QR-factorization  DF(xe)P  =  Q(R , 
ki  :=  0,  Vi  =  1,...  ,n; 

For  j  =  1, . . .  ,  to  Do:  If  ej"  =  P c'J*  Then  A:,  =  1; 

C  =  1; 

For  i  =  1 . n  Do:  If  ki  -  0  Then  zt  =  e?;  i  :=  t  +  1; 

Output  Z  :=  (zi,...  ,z»n)* 


3.  Sinipiicial  Approximations  and  Moving  Frames. 

Continuation  methods  are  probably  the  oldest  methods  for  the  computational  analysis 
of  an  implicitly  defined  manifold  M  although  this  fact  is  often  not  noticed.  They  apply 
to  problems  of  the  form  (1.1)  with  a  parameter  space  of  dimension  d  =  1  in  which  case 
the  solution  manifold  M  is  one-dimensional.  All  continuation  methods  begin  from  a  given 

point  x°  €  M  and  produce  a  sequence  of  points  x}.j  =  0. 1.2 . on  M.  In  general,  the 

step  from  x1  to  xJ+1  corresponds  to  an  implementation  of  a  local  parametrization  of  M. 
More  specifically,  a  one-dimensional  coordinate  subspace  T  =  span{u}  at  xk  is  chosen,  a 
predicted  point  w  is  selected,  and.  with  the  corresponding  local  coordinate  y  =  wTu  as 
input,  a  local  parametrization  algorithm,  such  as  GPHI  or  TPHI.  is  applied  to  generate 
the  next  point  x'  +  1  on  M.  The  various  continuation  methods  differ  in  the  choice  of  (i)  the 
coordinate  direction  u.  (ii)  the  construction  of  the  predicted  point  w.  and  (m)  the  form  of 
the  iterative  process  used  in  the  local  parameterization  algorithm. 

In  the  case  of  a  multidimensional,  implicitly  defined  manifold  A/,  it  is  obviously  diffi¬ 
cult  to  achieve  a  good  assessment  of  the  features  of  M  solely  from  information  along  some 
paths.  This  led  to  the  development  of  several  methods  for  the  approximation  of  subsets 
of  such  a  manifold  (see  jBrOGJ  for  a  recent  survey).  In  particular.  (RhSS)  presented  a  first 
method  for  computing  a  simplicial  approximation  S  of  a  neighborhood  of  a  given  point 

€  M  consisting  of  a  grid  of  points  xfc  €  A/.  t  =  0, 1 . *.  and  their  connectivity  pat¬ 

tern.  This  method  was  globalized  in  (BrRh94)  to  allow  for  the  computation  of  a  simplicial 
approximation  covering  a  specified  subset  of  M.  The  algorithms  in  [RHSS]  and  [BrRh94] 
were  restricted  to  the  case  d  =  2.  Recently,  in  [Br96],  an  extension  to  the  case  of  any  d  >  2 
was  developed  and  some  applications  and  other  related  algorithms  were  discussed. 

The  methods  in  the  cited  articles  are  similar  to  a  continuation  method.  At  a  point 
x>  on  the  current  “frontier''  of  the  already  computed  part  of  I  a  tangential  coordinate 
space  of  M  is  chosen  and.  by  means  of  a  tangential  local  parametrization  of  M  at  x;. 
the  already  computed  neighboring  simplices  of  I  incident  with  x'  are  mapped  onto  the 
tangent  space  Tz,  M.  In  addition,  a  reference  triangulation  patch  of  Rd  is  introduced  on 
Tz,  M  and  matched  to  the  open  facets  of  the  already  existing  part  of  H.  This  produces  a 
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set  of  points  in  the  unfilled  gap  of  the  triangulation  at  x>  which  are  then  used  to  complete 
the  triangulation  of  the  neighborhood  of  the  point  on  TtM.  While  for  d  =  2  this  is  a 
relatively  simple  task,  for  d  >  2  a  Delauney  triangulation  process  was  needed  in  [Br96]. 
The  nodes  of  the  completed  local  mesh  around  x1  are  then  mapped  onto  M  using  the  local 
parametrization  algorithm  TPHI  and  the  resulting  points  and  their  incidence  relations  are 

added  to  2. 

During  the  process  it  is  important  to  align  the  orientations  of  the  computed  simplices.  In 
the  setting  of  the  simplicial  approximation  algorithms  this  can  be  accomplished  by  setting 
the  orientation  of  any  newly  computed  simplex  equal  to  that  of  one  of  its  neighbors.  This 
was  the  approach  chosen  in  (Br96)  where  also  provisions  were  made  for  the  resolution  of 

any  conflicts. 

This  situation  is  a  special  case  of  the  more  general  problem  of  matching  the  basis 
orientations  of  the  local  parametrization  at  any  neighboring  points.  In  differential  geometry 
this  corresponds  to  the  concept  of  a  moving  frame;  that  is,  of  a  mapping  that  associates 
with  each'  point  x  of  a  d-dimensional  manifold  M  an  ordered  basis,  {ui, . . .  ,«<f}  of  TtM 

such  that  the  mappings  u,  =  1 . d,  form  d  vector  fields  on  M.  A  manifold 

is  parallelizable  if  such  a  moving  frame  exists  on  all  of  M. 

In  computational  problems  it  is  often  important  to  generate  a  moving  frame  on  an 
open  neighborhood  of  a  given  point  of  M.  Such  an  algorithm  was  developed  in  [Rh88]  and 
applied  in  the  mentioned  simplicial  approximation  method  for  two-dimensiional  manifolds. 

Under  Assumption  A,  suppose  that  orthonormal  bases  of  TSxM  and  TtiM  are  needed 
at  neighboring  points  xltx2  €  Af.  The  MANPAK  algorithm  COB  AS  uses  the  LQ- 
factorization  of  the  Jacobians  DF{x\)  and  DF(x2)  to  determine  the  desired  m  x  n  basis 
matrices  U\  and  U2,  respectively.  But,  as  noted  in  [CS84]  the  QR-factorization  algorithm 
(and  hence  the  LQ-factorization  algorithm)  need  not  produce  matrices  that  depend  con¬ 
tinuously  on  the  elements  of  the  given  matrix.  In  other  words,  the  computed  basis  U2  at 
xj  need  not  converge  to  the  basis  U\  at  xj  when  x2  tends  to  x\.  This  observation  extends 
to  other  algorithms  for  computing  the  nullspace  of  DF(x)  for  neighboring  x.  In  fact,  it 
relates  directly  to  the  well  known  loss  of  continuity  under  changes  of  the  matrix  elements 
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in  the  computation  of  eigenvectors  associated  with  a  multiple  eigenvalue. 

Let  A/0  be  an  open  subset  of  M  and  suppose  that  T  C  is  a  coordinate  subspace  at 
each  point  of  Mo-  Assume  that  Uq  is  an  orthonormal  basis  matrix  of  T  and  that  at  any 
x  €  Mo  an  orthonormal  basis  matrix  U(x)  of  TXM  is  chosen.  Let  V(x)  —  A(x)— B(x) 
be  the  singular  value  decomposition  of  the  d  x  H  matrix  V'(x)  =  U(x)rU0  and  form  the 
product  Q(x)  =  A(x)B(x)t  of  the  matrices  of  left  and  right  singular  vectors.  Then,  it  was 

proved  in  [RhSS]  that  the  mapping 

x  €  Mo  — *  U(x)Q(z)  €  £(RJ,Rn) 

is  of  class  Ck~l  on  M0  and  defines  an  orthonormal  moving  frame  on  M0. 

This  leads  to  the  following  MANPAK  algorithm  first  given  in  (R88): 

MOVFR  Input:  Orthonormal  tangent  basis  matrices  U0  and  Ut\ 

For  V  =  UjUo  compute  the  SVD  V  =  .4EBr; 

Q  =  ABt\ 

(Js  =  UtQ\ 

Output:  Rotated  basis  matrix  Ux. 

The  use  of  the  SVD  for  the  d  x  d  matrix  V  of  MOVFR  may,  of  course,  become  costly 
if  the  dimension  d  of  the  manifold  is  larger.  Thus  other  algorithms  for  matching  the 
orientation  of  computed  bases  are  of  interest.  A  simple  approach  is  based  on  the  use  of 
the  well  known  greedy  algorithm  (see  (PaRSS])  of  combinatorial  computing  which  can  be 
applied  to  two  orthonormal  bases  Ux  and  U2  at  some  neighboring  points  X|,  x2  6  M  in  the 

following  form: 

ORIENT  Input:  Orthonormal  basis  matrices  U\,  U2\ 

Compute  the  matrix  V  —  U^U2  =  (uij>  j  €  J ),  3  —  O’*** 
ui(i)  =0,  Vi  €  J-; 

For  i  =  1, . . .  ,  d  DO: 

Find  the  smallest  k  ~  k{  €  J  such  that  =  max{|u,;|,  j  €  J] 
and  u >(k)  =  0: 
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If  no  such  index  k,  exists  Then  re-orientation  failed; 
u{ki)  =  1; 

If  V{k  <  o  Then  change  the  sign  of  all  elements  in  the  k,- th  column  of  Ui\ 

End  If; 

With  the  permutation  P  €  Fef  —  t  *  €  3  form  U 2  :=  PU2 , 

Output:  Rotated  basis  matrix  1/ 2 . 

It  can  be  shown  that  the  algorithm  may  fail  in  cases  when  MOVFR  is  successful.  How¬ 
ever,  this  happens  very  rarely.  Moreover,  the  problem  can  be  resolved  by  incorporating  m 
ORIENT  a  standard  backtracking  approach  (loc.cit.).  This  modification  has  been  imple¬ 
mented  but  was  not  included  in  MANPAK. 

5.  Sensitivity  Computation. 

As  noted  in  the  Introduction,  many  physical  systems  lead  to  mathematical  models  in 
the  form  (1.1)  of  parameter  dependent  equations.  Then  it  is  often  assumed,  especially  in 
the  engineering  literature,  that  the  solutions  x  =  (z,  A)  of  (1.1)  can  be  written  in  the  form 
(z(A),  A).  Of  course,  under  Assumption  A,  this  requires  that  at  x  €  M  the  parameter 
subspace  A  is  a  local  coordinate  space  of  M,  which  means  that  ker  DF(x)  fl  A  =  {0}  or, 
equivalently,  that  the  partial  derivative  DzF(x)  =  DF{x)\z  of  F  with  respect  to  the  state 
space  is  nonsingular.  If  this  holds  then,  traditionally,  the  derivative 

(5.1)  Dz{  A)  =  -Dzfr(2(A),A)“1F((z(A),A) 

is  defined  as  the  sensitivity  measure  of  the  particular  solution  under  variation  of  the  pa¬ 
rameters  (see  e.g.  [TV79]).  As  noted,  (5.1)  is  applicable  only  at  solutions  x  €  M  which 
are  not  foldpoints  of  M  with  respect  to  A.  In  many  applications  this  is  an  undesirable 
restriction  since  exactly  these  foldpoints  are  of  special  interest.  In  fact,  it  is  one  of  the 
fundamental  observations  of  bifurcation  theory  that  these  are  the  points  where  the  char¬ 
acter  of  the  solutions  of  (1.1)  may  undergo  major  changes.  This  led  to  the  development  of 
a  more  general  sensitivity  theory  in  [Rh93]  which  applies  generally  on  M . 

Let  T  C  Rn  be  a  local  coordinate  space  of  M  at  a  given  point  x°  €  M  and  denote  by 
(Vd,  9)  the  induced  local  parametrization  of  M  at  xc.  Then  D<p( 0)y  represents  the  change 
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of  the  solution  xr  in  the  local  coordinate  direction  y  €  T.  Let  the  columns  of  the  n  x  m 
matrix  \'m  and  of  the  n  x  <1  matrix  V.{  form  orthonormai  bases  of  T  and  Tx,  respectively. 
Then  VmV*D?(to)y  is  the  orthogonal  projection  of  Dy>{ 0)y  onto  TL.  It  characterizes,  in 
essence,  the  change  due  to  the  nonlinear  nature  of  M:  in  fact,  if  M  were  Hat:  that  is.  if 
M  C  T,  then  this  vector  would  be  zero.  In  line  with  this,  in  j Rh93] .  the  linear  mapping 


(5.2)  £  €  £(Rd,Rn),  S  =  l'JlL-tO) 

is  defined  as  the  sensitivity  map  of  M  at  xc  with  respect  to  the  local  coordinate  space  T. 

In  [Rh93]  it  was  shown  that  when  the  natural  parameter  space  A  is  a  local  coordinate 
space  at  xc  €  A'/,  then  the  new  definition  (5.2)  reduces  exactly  to  (5.1).  More  generally, 
consider  besides  the  local  coordinate  space  T  at  xc  also  the  tangential  coordinate  space 
TZeM.  In  analogy  to  V'm  and  Vj,  let  the  columns  of  the  n  x  m  matnx  Um  and  of  the  n  x  d 
matrix  Ud  form  orthonormal  bases  of  Tte M  and  TZcM 1,  respectively.  Then  the  following 
relations  were  proved  in  [Rh03] 


(5.3)  UAv7u<r\ 

_  dist(r,Tr,A/) 

l5'4)  11-1,2  "  [1  -dist(T,rreA/)P  ' 

Here,  as  usual,  the  distance  between  any  two  equidimensional  linear  subspaces  5 \  and  5j 
of  Rn  is  defined  by  dist(S,,S2)  =  ||A  -  Pi ||2  where  P,.  Pi  are  the  orthogonal  projections 

onto  Si  and  S2,  respectively. 

In  MANPAK  two  algorithms  implement  (5.3)  and  (5.41.  respectively.  In  view  of  the 
noted  engineering  applications,  both  algorithms  assume  that  the  natural  parameter  space 
A  is  used  as  the  local  coordinate  space,  and,  more  specifically,  that  A  is  spanned  by  the 
canonical  basis  vectors  e"  ,  j  =  1,  •  •  •  ,  d  specified  by  a  given  index  set  J  C  {!,•••  > n}- 


SENMAP  Input:  Index  set  J  =  {*i,  -  - .  ,(<<},  orthonormal  basis  matrix  Ud  of  TZeM; 

Form  the  d  x  d  matrix  A  with  the  columns  Uj  e",  i}  €  J,  J  =  1 . 

Compute  the  LU-factorization  of  A: 

If  A  is  numerically  singular  then  Output:  ’’Undefined  sensitivity”. 
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Determine  the  indices  kj  £  J,  1  <  kj  <  n.  j  —  1, . . .  ,  m  —  n  —  d; 

For  j  =  1, . . .  ,  m  Do:  Solve  Awj  —  Uj ; 

Output:  S  :=  (uit, . . .  ,  u>m)T. 

SENNRM  Input:  Index  set  J  -  {*», . . .  ,  :</},  orthonormal  basis  matrix  Ud  of  TXcM; 
Form  the  d  x  d  matrix  A  with  the  rows  (e”  )rUd,  j  =  1, . . .  ,  d;; 

Compute  the  smallest  singular  value  crd  of  A ; 

If  Cd  =  0  Then  (  =  oo  Else  C  =  1; 

Output:  ||E||2  :=  C- 

The  corresponding  algorithms  for  the  more  general  case  of  an  arbitrary  local  coordinate 
space  T  were  not  included  in  MANPAK  due  to  their  infrequent  applicability. 

6.  Derivatives  of  Local  Parametrizations. 

Let  M  be  a  d-dimensional  Ck  submanifold  of  RB  with  k  >  2  and  a  local  Ck 

parametrization  of  M  near  xc  €  M,  Then,  as  noted  in  Section  2,  for  any  ve  €  R  such 
that  (xc,ue)  €  TAf,  the  pair  (V*'  x  Rd,  (v p,D<?))  is  a  C*”1  local  parametrization  of  the 
tangent  bundel  near  { xe,ve ).  This  means  that  for  the  evaluation  of  the  corresponding 
local  parametrization  of  TM  at  the  that  point  we  need  an  algorithm  for  computing  the 
derivative  Dy  of 

Under  Assumption  A,  let  T  be  a  coordinate  subspace  of  M  at  xc.  As  before,  suppose 
that  on  Rn  and  R*  the  canonical  bases  are  chosen  and  that  the  columns  of  the  n  x  m  matrix 
U  form  an  orthonormal  basis  of  the  given  coordinate  subspace  T  at  xc.  Then  by  Theorem 
3.1  the  mapping  (3.2)  is  a  local  diffeomorphism  near  xc  6  M  and  the  local  parametrization 
(V*,v>)  induced  by  T  satisfies  (3.3).  Thus,  by  differentiation  of  (3.3),  it  follows  that 

(6.1)  DH(v(y))Dv{y)v  =  Jv ,  Vy  €  V*,  v  €  RJ. 

Since  the  Jacobian  (3.5)  of  H  is  nonsingular  at  xey  this  shows  that  at  any  x  =  ?(y),  V  €  Vd , 
the  derivative  Dip(y)  cam  be  computed  as  follows: 
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DGPHI  Input:  Local  coordinate  basis  U  at  Jacobian  DF(x)  at  a  neighboring  point 

x  =  r(y); 

Compute  the  LU-factorization  of  A  := 

For  j  =  1. .  ■  •  ,  d  Do:  Solve  Az >  =  e'j+J; 

Output:  D<?(y)  :=  (zj, . . .  , zj). 

In  line  with  our  general  format,  the  algorithm  assumes  that  the  Jacobian  DF(x)  has 
already  been  evaluated  by  the  calling  program. 

In  the  case  of  a  tangential  coordinate  space  we  obtained  the  modified  algorithm  TPHI 
based  on  the  LQ-factorization  of  DF(xe).  This  reduced  the  cost  of  the  GPHI  algorithm. 
For  the  evaluation  of  the  Jacobian  of  *  the  use  of  the  LQ-factorization  is  not  as  advan¬ 
tageous.  Suppose  again  that  we  want  to  compute  D<p{y)  at  x  =  ?(v)-  Moreover,  let  Ux 
and  Utc  ben  xd  matrices  with  orthogonal  columns  that  form  bases  of  TXM  and  TXtM, 
respectively.  With  the  matrix  representation  of  DH(x)  used  in  DGPHI,  the  j-th  column 
Zj  of  Dtp(y)  satisfies  DF(x)zj  =  0  and  UXtZj  =  e<j.  Hence  we  have  zj  =  Uxyj  for  some 
yj  €  Rd  and  UjUxyj  -  tj  which  implies  that 

(6.2)  D^y)  =  Ux[UjUx)-'. 

In  analogy  to  DGPHI  this  gives  the  following  algorithm: 

DTPHI  Input:  Tangent  basis  UXt  at  xc,  a  neighboring  point  x  €  M; 

Use  COBAS  to  evaluate  Ux  at  x; 

A  :=  U£U'i 

For  j  =  1, . . .  ,  d  Do:  Solve  Az,  =  ej;  z}  :=  Uxz}\ 

Output:  Dy>{y)  :=  (zi ,  •  •  •  <  =d)- 

Since  the  computation  of  the  basis  matrix  U,  requires  the  application  of  COBAS,  the 
use  of  DTPHI  is,  in  general,  more  costly  than  that  of  DGPHI. 

Suppose  now  that  k>  2  in  Assumption  A.  Then  9  is  at  least  tw.ce  differentiable  at  any 
t/  €  Vd  and  by  differentiation  of  (6.1)  we  obtain 

(6.3)  DH(<f{y))D2v{y)(vx,v2)  =  -D2  HMy^D^vx^D^vi ),  €  R  • 


DF(x)\ 

u t  ;; 
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which,  because  of  the  nonsingularity  of  DH(<f(y)),  defines  D2v(y)( uniquely.  In 
line  with  our  reverse  communication  paradigm,  the  following  MANPAK  algorithm  of  D2<p 
assumes  that  the  vector  w  =  D2F(x){uuu7)  has  already  been  computed  for  given  x  =  V»(y) 
and  u,  =  Dtf(y)vii  i  =  1,2. 


D2GPHI  Input:  Tangent  basis  Ut ,  at  xc,  Jacobian  DF(x)  at  a  neighboring  point  x  €  M. 
the  vector  w  =  D2F(x)(ui,U2 )  with  Uj  =  D<p(y)vi,  i  —  1,2; 

Compute  the  LU-factorization  of  A  I  ^  1  i 

/  -u>\ 

Solve  Az  =  for  z; 

V  o  ) 

Output:  D2(p(y){v\,V2) z. 


7.  Curvature  and  the  Second  Fundamental  Tensor. 

A  principal  application  of  the  algorithms  for  the  computation  of  the  derivatives  of  a 
local  parametrization  arises  in  connection  with  the  solution  of  differential  algebraic  equa¬ 
tions.  This  is  discussed  in  the  mentioned  companion  paper  [Rh96c]  and  will  not  be  further 
addressed  here. 

Another  interesting  use  of  the  algorithm  DGPHI  occurs  in  certain  constrained  minimiza¬ 
tion  methods.  We  indicate  here  only  briefly  the  approach  suggested  in  [Rh96a].  Under 
Assumption  A  with  k  >  2,  suppose  that  j  :  £  6  R"  h  R1  is  a  Cr  functional,  r  >  2, 
on  some  open  set  E  and  that  E  n  M  is  not  empty.  Consider  the  problem  of  computing 
a  local  minimizer  x*  €  M  of  g  on  M,  and  let  x  €  M  be  a  point  on  M  that  represents 
our  current  approximation  of  x\  We  introduce  a  local  parametrization  (V*,?)  of  M  at 
x.  Then  the  local  representation  h  =  goyoig  near  x  is  at  least  of  class  C 2  and,  locally, 
the  minimization  of  g  on  M  is  equivalent  with  the  unconstrained  minimization  of  h.  This 
suggests  the  application  of  a  trust  region  step  to  h  in  order  to  obtain  a  new  approximation 
of  x*.  For  this  we  approximate  h  by  the  quadratic  functional 


(7.1) 


hq(y)  =  h(0)  +  Dh( 0)y  +  ~ D2h(0)(y,y ),  y  6  Rd. 
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Here  we  have 


OHO)!,  =  09(*)D9( 0)»,  D‘hmv,v)  =  0j(I)0,9(O)(y,9)  +  D'gWDmv.D'Xm. 


which  can  be  evaluated  by  means  of  the  MANPAK  routines  DGPHI  and  D2GPHI.  For 
details  about  the  trust  region  step  we  refer,  for  example,  to  |MS81|.  Clearly,  in  practice, 
the  computation  of  the  local  Hessian  Dh( 0)  can  be  replaced  by  some  update  scheme. 


The  second  derivative  D\  of  a  local  parametrization  has  an  important  connection  with 
the  second  fundamental  tensor  of  the  manifold  U.  This  tensor  is  a  concept  of  Riemanman 
geometry;  that  is,  it  re.,uires  a  metric  on  the  manifold  M.  We  use  here  the  metric  induced 
by  the  canonical  (Euclidean)  inner  product  of  R".  For  a  definition  of  this  symmetric. 

vector  valued  tensor 


(7.2)  V, :  T,M  x  T,M"  T.M1,  z€:V, 

in  a  setting  similar  to  that  used  here,  we  refer  to  (RaRh90|.  In  lieu  of  a  definition  we  cite 
only  the  following  characterization  of  Vt  proved  in  [RaRh95j. 

Theorem  T.l.  Under  Assumption  A,  let  i  €  M,  Z  any  complement  of  T.M,  and  Q 
the  orthogonal  projection  onto  (T.M)1.  Then  the  component  of  the  second  fundamental 
tensor  V,  of  M  at  r  in  the  tangential  directions  u‘,u3  S  T.M  ts 

,7.21  P,(u1,uJ)--<3(Df'(s)|zr,0Jf(*X“‘.“i)- 


Let  T  be  a  local  coordinate  space  at  i  €  M  and  ( V,  v>)  the  induced  local  parametrization 
of  M.  Moreover,  set  Z  =  T1,  and,  as  before,  let  the  columns  of  the  n  x  d  matrix  U  form 
an  orthonormal  basis  of  T.  Note  that  at  z  =  <p(0)  the  equation  (6.3)  defining  DM 0) 
requires  that  DF(.)DV0)(»h»i)  -  -D’Wfu’.u3)  and  U1-DV(0)(v„«)  =  0  with 
u,-  =  D<?{0)vi,  i  =  1,2.  Evidently,  this  is  equivalent  with 

[DF(x)\Z]D\(0)(vi,v2)  =  - D2F(x)(u'.u 2), 


whence,  it  follows  from  (7.2)  that 

V’.(u,.u2)  =  QZ>V(0)(v,,v2),  =  Dv(0)vi,  i  =  1,2. 


(7.3) 
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Note  that  £>V(0)(vi, ”2)  €  T1,  and  thus  V*(u\u2)  =  £2y5(0)(ui' ^2)  for  the  tangential 
coordinate  space  T  =  TtM.  Hence  this  tensor  component  can  be  obtained  by  application  of 
the  MANPAK  algorithm  D2GPHI  with  a  tangential  coordinate  space.  However,  if  COBAS 
was  used  with  ZT  =  DF{ xc)  to  compute  the  tangent  basis,  we  may  proceed  as  in  DTPHI 

and  simplify  the  process  as  follows: 

TSFT  Input:  The  LQ-factorization  DF(x)  =  PT{L  0)QT, 

the  vector  w  =  I?2.F(x)(ui,  tij)  with  u<  =  D<p(y)vi,  i  —  1,2, 

Solve  Lz  =  PTw  for  z  €  Rm; 
z  :=  Q\z\ 

Output:  V*,(vi,t>2)  :=  z- 

The  second  fundamental  tensor  characterizes  curvature  properties  of  the  manifold  and 
is  closely  connected  with  the  Riemann  curvature  tensor  R  of  M.  Accordingly,  it  is  not 
surprising  that  it  has  numerous  applications.  For  example,  in  [RaRh90]  the  tensor  has 
been  applied  for  computating  bifurcation  directions  at  certain  foldpoints  on  an  implicitly 
defined  manifold.  Then,  in  [RaRhD5]  it  was  used  in  an  algorithm  for  solving  the  Euler- 
Lagrange  equations  arising  in  the  modeling  of  constrained  dynamical  systems.  In  both 
cases,  early  forms  of  the  MANPAK  algorithms  D2GPHI  and  TSFT  were  applied. 

Note  that  it  suffices  to  have  a  method  for  computing  the  diagonal  components  Vx(u,u) 
since,  by  the  bilinearity  and  symmetry  of  the  tensor,  we  have 

Vt(ul,u2)  =  2Ft(u,u)-i(Vt(u1,u1)  +  Vl(u2,u2)) 

for  any  u1 ,  u2  €  Tr,  M  and  u  =  (l/2)(u1  +  u2).  In  (RR90a]  a  geometrically  based  algorithm 
was  given  for  approximating  the  diagonal  component  Vr(u,u)  for  any  u  €  TtM. 

Let  it  be  any  path  on  M  and  xt,xe,xr  €  M  three  consecutive  points  along  tt  that  form 

a  nondegenerate  triangle.  Then  the  curvature  of  the  circumscribing  circle  of  the  triangle 

is  given  by  Heron’s  formula 

4  / _  — - - ^  a  +  b  +  c 

(7.4)  k  =  —  y/s(s  -  a){s  -  b)(s  -  c),  s  =  - 
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where  a  =  ||x,  -  xr||2,  b  =  ||xr  -  x,j|„  c  =  ||x,  -  x,||r  When  the  outside  points  x,,xr 
tend  to  the  middle  point  x,  then  the  circumscribing  circle  tends  to  the  osculating  circle  of 
the  path  rr  at  xr  and  the  limit  of  k  is  the  curvature  of  tr  at  xe.  It  was  shown  in  [RR90a] 
that  the  value  (7.4)  of  k  approximates  ||V^(u.u)||2  for  the  tangent  vector  u  €  Ttt M  of 
-  at  Moreover,  the  unit  vector  in  the  direction  of  V'rt(u.u)  is  approximated  by  the 

principal  normal  vector  of  the  path  ~  at  xr. 

For  the  computation  it  is  useful  to  rewrite  (7.4)  as 

,  ,  ,  , _  1 

k  =  -  ( - hr]  V 1  -  <!>2  sina,  6  =  ae  -  be,  a  -  arcos  7,  7  -  .  ^  • 

C  t?c 

This  leads  to  the  following  MANPAK  algorithm: 

CURVT  Input:  Consecutive  points  x<,xc,xr  along  a  path  tt  on  M. 

tangent  vector  v  €  TZrM  of  tt  at  xc,  machine  precision  e 
Evaluate  a  :=  ||x<  -  xc||2,  b  :=  j|xr  -  xd|2,  c  :=  ||xf  —  xr||2, 

“  :=  v/IMI*! 

p  :=  ( 1/ a)(i  -  x,)  +(l/6)(xr  -  Xc); 
p  :=  p  -  (pTu)u; 
ar  :=  a/c;  bc  :=  b/c\ 
b  :=  ae  -  bc\  7  :=  1/1"  h  br)\ 

If  1  -  7  <  <  Then 

k  :=  0; 

Else 

a  :=  arc  cos  7; 

k  :=  (l/c)(l/ac  +  1/M  v'l  -  sina; 

End  If; 

Output:  p,  k;  KIc(u,u)  :=  *p. 

When  1  -  7  falls  below  the  machine  precision  «  then,  in  floating  point  arithmetic,  a  will 
be  zero  and,  accordingly,  we  set  k  =  0.  The  approximation  p  of  the  principal  normal  of 
the  path  at  xc  is  generated  by  a  simple  Gram-Schmidt  orthonormalization  step. 

For  some  numerical  examples  involving  CURVT  we  refer  to  (RROOa). 
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